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We compute the gravitational self-force (or "radiation reaction" force) acting on a particle falling 
radially into a Schwarzschild black hole. Our calculation is based on the "mode-sum" method, in 
which one first calculates the individual ^-multipole contributions to the self-force (by numerically 
integrating the decoupled perturbation equations) and then regularizes the sum over modes by 
applying a certain analytic procedure. We demonstrate the equivalence of this method with the 
(^—function scheme. The convergence rate of the mode-sum series is considerably improved here 
(thus reducing computational requirements) by employing an analytic approximation at large I. 
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The space-based gravitational wave detector LISA 
(Laser Interferometer Space Antenna), scheduled for 
launch around 2011 jjj, will open up a window for the low 
frequency band below 1Hz, allowing access to a variety 
of black hole sources. A main target for LISA would be 
the outburst of gravitational radiation emitted during the 
capture of a compact star by a supermassive black hole — 
a 10 5-7 solar masses black hole of the kind now believed 
to reside in the cores of many galaxies, including our own 
0j . By LISA's launch time, a sufficient theoretical under- 
standing of the orbital evolution of such systems, includ- 
ing radiation reaction effects, must be at hand, to allow 
design of accurate templates necessary for detection and 
interpretation of the gravitational waveforms. Due to the 
extreme mass-ratio typical to the binary systems of in- 
terest, the entire problem can be conveniently treated in 
the context of perturbation theory, being a relatively ma- 
ture branch of gravitational physics: the compact stellar 
object is thus modeled by a point-like particle, and its 
field treated as a perturbation over the fixed Kerr geom- 
etry of the large black hole. To leading order in the mass 
ratio, such a particle then traces a geodesic of the back- 
ground spacetime, and one asks about radiation-reaction 
induced corrections to this geodesic. 

The "traditional" approach to this problem relies on 
energy-momentum balance considerations (|]: By calcu- 
lating the fluxes of energy and angular momentum to in- 
finity and across the event horizon, one attempts to infer 
the temporal rate of change of the particle's "constants" 
of motion. This technique is applicable only in adiabatic 
scenarios, in which the time scale for radiation reaction 
effect is much larger than the dynamic time-scale of the 
system; it is not clear yet whether this approximation 
is valid for the entire range of relevant LISA parame- 
ters Q. Moreover, this approach seems insufficient for 
tackling generic orbits in Kerr spacetime (i.e., ones both 
eccentric and inclined) even under the adiabatic approx- 
imation. This has led many researches, particularly over 
the last five years, to turn to the useful notion of the local 



self-force (SF). 

Consider a point-like particle of mass /i moving freely 
around a black hole with mass M 3> /z; and treat the 
particle's gravitational field, h a p oc 0(/x), as a linear per- 
turbation over the background metric g a p. In the "SF 
picture" , this particle's equation of motion is written as 

/i^V 3 = F£f , (1) 

where u a is the particle's four-velocity, a semicolon de- 
notes covariant differentiation with respect to g a p, and 
F s " lf oc 0(a* 2 ) describes the leading-order SF effect. Since 
the perturbation h a p obviously diverges at the particle's 
location, the problem of obtaining involves the in- 
troduction of a reliable regularization scheme. A well 
established formal prescription for constructing F^ lf , re- 
lying on a physically consistent regularization method, 
became available recently with the work of Mino, Sasaki 
and Tanaka (MST) pj. The same formal prescription 
was introduced independently by Quinn and Wald (QW) 
II, based on an axiomatic approach. To allow a practi- 
cal implementation of the MST/QW formal prescription 
in actual calculations, a "mode-sum" scheme was later 
devised in Refs. §, Q, §, based on the MST/QW re- 
sult. An alternative regularization approach, based on 
the £— function technique, was introduced in Ref. ||. 

The main objective of this paper is to report on a 
first actual calculation of the gravitational SF, based 
on the mode-sum prescription. Focusing, as a test-case, 
on radial trajectories in a Schwarzschild background, we 
demonstrate the applicability of this approach, and push 
forward some analytic and numerical techniques which 
may later be applied to more general orbits. Among 
the new results presented here: (i) Two different deriva- 
tions of the "regularization parameters" , independent 
of each other and of the derivation of ||; (ii) Consis- 
tency of the MST/QW regularization with the ^-function 
method; (iii) A first explicit example of the gauge- 
invariance feature of the regularization parameters pre- 
dicted in Ref. (Hi; (iv) An improved numerical method 
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for integrating the decoupled field equations to fourth or- 
der accuracy; (v) An analytic approximation developed 
for improving the convergence rate of the mode-sum se- 
ries (see below). Full details of our analysis shall be pro- 
vided elsewhere pL 

Throughout this paper we use "geometrized" units 
G = c = 1, metric signature — H — h, and the standard 
Schwarzschild coordinates t, r, 9, ip. We consider a parti- 
cle falling radially into a Schwarzschild black hole with 
mass M, starting at rest at r — tq, and let r p denote the 
value of r at the SF evaluation point. In the lack of SF, 
the particle traces a geodesic characterized by the (con- 
served) specific energy E = (1 — 2M /ro) 1 / 2 . By virtue 
of the symmetry of the above setup, we obviously have 



F p lf = 0; we hereafter thus focus on the r, t com- 



ponents of the SF. 

Let us start by briefly reviewing the mode-sum method 
for constructing the gravitational SF: First, one has to 
calculate the multipole modes of the metric perturba- 
tion h a p in the harmonic gauge, denoted here by h aB 
(this refers to the quantity obtained by summing over all 
azimuthal numbers m and over all ten tensor harmon- 
ics for a given multipole number £). This calculation 
is done through a numerical integration of the decou- 
pled linearized Einstein equations. Then, one constructs 
the ^-mode contribution to the "full" force, denoted here 
by F ai , through a certain operation involving lst-order 
derivatives of h l a g [see Eq. (15) of Ref. jsj]. In the radial 
motion case, this operation reduces to 



pat 
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(evaluated at the particle's location), where h a p = 
hapg^h' k a ^ s = u^g aS /2 + g^g a5 /A + 



"a/3 

u a g^u s /4:-g a f 3 u'u s -u a uf 3 u^u s /2, and the ± sign cor- 
responds to taking the derivative from r — > r p , respec- 
tively. (Note that these force-modes satisfy the normal- 
ization condition u a F± l = 0.) While the perturbation 
itself diverges at the particle's location, the individual 
modes h^p are continuous everywhere || — an important 
benefit of the mode-sum approach (this holds in the har- 
monic gauge or in any other gauge related to it by a 
regular gauge transformation). Typically, however, the 
derivatives of h a p are found to have a finite discontinu- 



ity through the particle's location, yielding two different 
finite values F± e . According to the mode-sum method, 
the gravitational SF is then constructed through H 
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A a ± L 



B a _ c a /L] - D a , (3) 



where L = £ + 1/2 and the (^-independent) quantities 
A a , B a , C a , and D a are the so-called "regularization 
parameters", whose values depend on the orbit under 
consideration. Roughly speaking, the expression A±L + 
B a + C a /L reflects the asymptotic form of F± e at large 



£ [ensuring convergence of the sum in Eq. while the 
parameter D a is a certain residual quantity that arises in 
the summation over £. (See 0, |^[ for an exact definition 
of these parameters.) 

Incorporating a systematic perturbation expansion of 
the £— mode Green's function associated with the pertur- 
bation equations in the harmonic gauge — an implementa- 
tion of the technique developed in [Q — we have obtained 
for radial trajectories JlT|, 
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Z3*> B=-^^> (4b) 



C a = D a = 0, 



(4c) 



where / = 1 - 2M/r p and r p = -(E 2 - f) 1 ' 2 . These 
values agree with those derived (for generic orbits) in 
Ref. H using a different method. Note that whereas the 
values of A±, B a , and C a can be verified numerically by 
examining the behavior of the modes F± at large £ (see 
below), the value of D a cannot be so verified; hence the 
importance of our independent derivation of D a . 

The above prescription requires one to tackle the per- 
turbation equations in the harmonic gauge. These equa- 
tions are separable with respect to £, m , but it is not 
clear how, or whether at all, one could avoid, the cou- 
pling occurring between different elements of the tensor- 
harmonic basis. A more practical derivation of h a g (and, 
consequently, of F± e ) is possible in the Regge- Wheeler 
(RW) gauge (l^, 0|: Here, a well developed formal- 
ism [D| allows one to derive all h aQ components from 
two scalar generating functions, by mere differentiation 
H O. These two waveforms satisfy a scalar-like wave 
equation which is easily accessible to numerical treat- 
ment. Now, it has been shown jl(| that the mode-sum 
formula ^ is valid, with the same parameter values, for 
any gauge related to the harmonic gauge through a regu- 
lar gauge transformation. The RW gauge belongs to this 
regular family of gauges so long as radial trajectories are 
considered [l0[ | — as in our current work. This shall allow 
us to work here entirely within the convenient RW gauge. 

Using a variant of the Green's function expansion tech- 
nique mentioned above — this time applied to the pertur- 
bation equations in the RW gauge — we have been able 
to directly obtain the values of A a , B a , and C a associ- 
ated with the RW-gauge modes F^° . (The details of this 
derivation shall be given in 0.) The RW-gauge param- 
eters thus obtained were found to have, in the head-on 
case considered here, precisely the same values as in the 
harmonic gauge. This serves as a first explicit demonstra- 
tion of the regularization parameters' gauge-invariance 
property predicted in p0[ . 
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It should be noted that the mode-sum prescription (j|) , 
stemming from the standard MST/QW regularization 
scheme, completely conforms with the ^-function regu- 
larization approach introduced in Ref. In the latter 
too, the SF is brought to the form (||), with the param- 
eter D a shown ||, Q to be oc ((0, 1/2) = (where ( is 
Riemann's generalized zeta function). 

Incorporating the parameter values (Q) in Eq. (||) 
(noticing that A+ = —A"), we next write the mode-sum 
formula in the compact form 



0.0 F 



E 

£=0 



(5) 



where F al = (Ff + F^)/2, and F r f g [(x 0(£- 2 )} are 
the "regularized" modes. Note that B a {r p ) describes the 
asymptotic form of F al at £ — > oo. 

We now turn to the actual implementation of the pre- 
scription (|J), beginning with the numerical calculation 
of F al . As already mentioned, our calculation was car- 
ried out within the RW gauge. All tensorial components 
of the €-mode metric perturbations in the RW gauge (in 
fact, only the even-parity part of h a p plays role in our 
head-on case) are conveniently constructed from a single 
scalar generating function — Moncrief's gauge-invariant 
waveform ip* . This construction, prescribed in || , in- 
volves twice differentiating i//. Then, the desired modes 
F ag are obtained using Eq. (0). Thus, the numerical task 
reduces to integrating the (inhomogeneous) hyperbolic 
wave equation satisfied by ^ l fl6fl , with the appropriate 
source term associated with the point-like particle, and 
with a proper choice of initial data. This numerical prob- 
lem has been formalized and worked out previously in 
Rcf . [[l6| . We have developed an improved version of the 
above numerical scheme, which ensures fourth-order nu- 
merical convergence. (This is essential for our purpose, 
as the construction of F al involves three derivatives of 
the numerical- integration variable ip .) 

Typical results from applying the above numerical pre- 
scription are presented in Fig. [l], showing the first few 
(averaged) modes F as a function of r p for the sample 
value r = 14M (E = 0.926), and demonstrating the an- 
ticipated oc £~ 2 behavior of F^ g . The above construction 
of F is only applicable to the "radiative" modes I > 2. 
The modes I = 0,1, although merely reflecting a residual 
gauge freedom 113] , must also be taken into account in 
the mode sum ^). As shown in ]T^ ], the I = 1 even- 
parity perturbation is completely removable by a gauge 
transformation (interpreted as a translation to the cen- 
ter of mass system), and we may take F^ = 0. The 
I = perturbation mode (interpreted as a variation in 
the mass M) is constructed analytically in jllj] — the re- 
sultant contribution F^° is also plotted in Fig. |l|. [jlTj 

Since F r "| oc 0(£~ 2 ), the mode sum in (|J) admits the 
slow convergence rate ccl/£. This means that achieving 
even a modest accuracy requires one to sum over many 
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FIG. 1: Upper figure: The "full" modes F rt {r v ) at I = 
0, . . . ,8, for a particle released from rest at ro = 14M. 
Note how a limiting curve [given analytically by B r (r)] is 
approached at large I. The wavy feature near r p = 14M is 
due to the radiation content of the initial data (chosen here 
as conformally flat). This "spureous" feature damps down by 
the time the particle reaches r p ~ 10M, exposing the inherent 
self force effect. The bottom figure demonstrates the antici- 
pated oc £~ 2 convergence of the difference Ff* s = F rl — B r . 



modes, which is numerically very demanding (numerical 
integration of the decoupled field equations becomes in- 
creasingly difficult at growing £) . To improve the conver- 
gence rate of the mode-sum series, we have obtained O] 
an analytic approximation for F£~^ at large £: By extend- 
ing the local analysis of the £— mode Green's function one 
order beyond the calculation of the three parameters A a , 
B a , and C a , we have obtained an analytic expression for 
the 0(L~ 2 ) term of F ai . This significantly improves the 
mode-sum convergence, especially by virtue of the fact 
that the next, 0(L~ 3 ) term in the mode-sum is expected 
to vanish (this can be shown for any positive odd power 
of 1/L in the mode sum using straightforward parity ar- 
guments H], and is further supported by our numerical 
results — see Fig. 0). We have found [O 
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where r is the particle's proper time (the two expressions 
are negative definite for all r < ro). These expressions 
are in perfect agreement with the numerical results, as 
demonstrated in Fig. |^. The contribution of the 0(L~ 2 ) 
expansion term to the overall SF is now easily obtained 
analytically, using J2h=2 ^ 2 — 0-49. The remainder of 
the mode sum now converges as oc £~ 3 . By calculating 
numerically only the first 10 modes (say), one now ob- 
tains the SF to within a mere relative error of ~ 1CP 3 
(compare this with a ~ 10 _1 error when not using the 
analytic approximation) . 
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FIG. 2: Analytic approximation vs. numerical results: The 
plot shows the numerically calculated (regularized) modes 
-Frog along with their large-^ analytic approximation, Fan a i yt j c , 
as given by the 0(L~ 2 ) term in Eq. (jfia]). Shown are the 
modes £ — 0, 1, 2, 3, 4 for ro = 1AM [times L 2 for the sake of 
comparison] . The inset shows the reminder F^g — F ar f alytic at 
r p = 6A/, demonstrating its anticipated oc £~ 4 behavior. The 
wavy feature at the onset of the plunge is associated with the 
"spurious" radiation content of the initial data; the inherent 
SF is exposed only after these waves are dissipated away. 



Figure || shows both r and t components of the overall 
SF resulting from summing up all individual mode contri- 
butions. The modes £ — 2, . . . ,8 were obtained numeri- 
cally, while for I > 8 we used the analytic approximation 
(0) (for £ = 0, 1 we used the exact solutions mentioned 
above). The radial component of the SF is found to point 
inward (i.e., toward the black hole) throughout the entire 
plunge. This seems to be a universal feature which does 
not depend on the starting point tq. Consequently, the 
work done by the SF on the particle is positive, resulting 
in that the energy parameter E increases throughout the 
plunge. The instantaneous rate of change of E is given 
by @ fJ-(dE/dr) = -F t (> 0), and the total change of 
E is obtained by integrating this expression along the 
worldline from r(ro) to r(r = 2M). It is important to 
stress, however, that this result will be attached to our 
specific choice of gauge (as opposed to the energy flux at 
infinity, which is gauge invariant) Jl9| . 

In summary, the mode-sum approach for calculating 
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FIG. 3: The upper and bottom panels show, respectively, 
the r and t components of the overall SF on a particle start- 
ing at rest at ro = 14M. The plots labeled as "numerical" 
are produced by summing up the numerically calculated (reg- 
ularized) modes up to £ = 8, and then incorporating our an- 
alytic approximation at £ > 8 (these higher modes contribute 
up to ~ 20% of the total force). Also given, for compari- 
son, is a curve based entirely on the analytic approximation 
(0) (summed over £ = 2, . . . ,oo plus the exact solutions for 
I — 0, 1); and a curve showing the mere 1 = 2 contribution. 
The latter serves to illustrate the importance of higher £ con- 
tributions. All curves reach a finite value at the horizon. 



the gravitational SF was successfully applied here in 
the test case of radial motion in Schwarzschild space- 
time. We have also demonstrated the feasibility of apply- 
ing an analytic approximation for improving the mode- 
sum convergence and even providing a rough estimate 
to the SF. This marks a significant milestone in our 
(still long) way toward being able to compute the or- 
bital evolution of generic orbits in Kerr spacetime. The 
next step along this way, already being considered, is the 
implementation of the mode-sum prescription to more 
general orbits in Schwarzschild background. This will 
provide a first opportunity for validating the SF ap- 
proach against available calculations based on the stan- 
dard energy-momentum balance approach. 
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